Abstract. We have combined time dependent hydrody- 
namics with a two-fluid model for dust driven AGB winds. 
Our calculations include self-consistent gas chemistry, 
grain formation and growth, and a new implementation of 
the viscous momentum transfer between grains and gas. 
This allows us to perform calculations in which no assump- 
tions about the completeness of momentum coupling are 
made. We derive new expressions to treat time dependent 
and non-equilibrium drift in a hydro code. Using a sta- 
tionary state calculation for IRC -f 10216 as initial model, 
the time dependent integration leads to a quasi-periodic 
mass loss in the case where dust drift is taken into ac- 
count. The time scale of the variation is of the order of 
a few hundred years, which corresponds to the time scale 
needed to explain the shell structure of the envelope of 
^ ; IRC -hl0216 and other AGB and post-AGB stars, which 
. has been a puzzle since its discovery. No such periodicity 
P\J ' is observed in comparison models without drift between 
dust and gas. 
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1. Introduction 



Dust dit'lvoii wliidH are powered by a raHcliiaLliig liiLei'play 



tion of grain nuclcation and growth, have been developed 
in the Berl i n group, initially for ca rbon-rich objects ( Gail 



of radiaLiuii, diemical I'uacLiuiiy, sLellar pulsaLiuiis and dy- 
namics. As soon as the envelope of a star on the Asymp- 
totic Giant Branch (AGB) develops sites suitable for the 
formation of solid "dust" (i.e. sites with a relatively high 



et al. 1984 ; Gail fc Sedlmayr 1987 ) and more recently also 



for the more com plicated case of silicates in circumstellar 

shells of M stars (|Gail fc Sedlmayr 19991) . 

Two-fluid models, in which dust and gas are not neces- 



densitji and a low LempeiaLuie) iLs dynamics will be dum- 
inated 



by ladiaLiuu piessuie. DusL giaiiis aie exLiemely — ^ 
sensiti'je Lu Llie sLellai ladiaLiuu and expeiieuce a laige 
radiation pressure. The acquired momentum is partially 
transferred to the ambient gas by frequent collisions. The 
gas is then blown outward in a dense, slow wind that can 
reach high mass loss rates. 

The detailed observations of (post) AGB objects and Plan- 
etary Nebulae (PN) that have become available during the 
last decade have shown that winds from late type stars are 
far from being smooth. The shell structures found around 



(1992 



sarily co-moving, have been less well studied. Berruyer & 
Frisch (19831 , jBerruyer (1991)| and [MacGregor fc Sten- 



pointed out that, for stationary and isother- 
mal envelopes, the assumption of complete momentum 
coupling breaks down at large distances above the pho- 
tosphere and for small grains. Self-consistent, but again 
stationary, two-fluid models, considering the grain size 
distribution, dust formation and the radiation field were 



developed by Kriiger and co-workers (Kriiger et al. 1994 



Kriiger fc Sedlmayr 1997) 



e.g. CRL 2688 (the "Egg Nebula", Ney et al. 1975; Bahai 



et al. 199g), NGC 6543 (the Cat's Eye Nebula, Harrington 



The only studies in which time dependent hydrodynam- 
ics and two-fluid flow have been combined so far are the 
work of 



Mastrodemos et al. (1996) 



fc Borl t.uwski 1994) and Lliu AGO sLai IRQ -H0216 (Me 



ron fc Iluggins 1999, 2000), indicaLe LhaL Lhe uuLfluw has 
quasi-periodic oscillations. The time scale for these oscil- 
lations is typically a few hundred years, i.e. too long to be 
a result of stellar pulsation, which has a period of a few 
hundred days, and too short to be due to nuclear ther- 
mal pulses, which occur once in ten thousand to hundred 
thousand years. 

Stationary models, in which gas and dust move outward 
as a single fluid, do not suffice to explain the observa- 
tions. Instead, time dependent two-fluid hydrodynamics, 
preferably including (grain) chemistry and radiative trans- 
fer, may help to explain the origin of these circumstellar 
structures. 

Time dependent hydrodynamics has been used t o study 
the influence of stellar pulsations on the outflow ( Bowen 



Schonberner 2000) 



and that of the Pots- 



dam group (^teffen et al. 1997| ; [Steffen et al. 1998| ; [Steffen 



198 j ; Fleischer et al. 1992). The coupled system of radi- 



ation hydrodynamics and time dependent dust formation 
was solved by Hofner et al. (1995)| . 

Stationary calculations, focused on a realistic implementa- 
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In the next section, we will argue that time dependence 
and two-fluid flow are not just two interesting aspects of 
stellar outflow but that they have to be combined. It turns 
out that fully free two-fluid flow, i.e. in which no assump- 
tions at all about the amount of momentum transfer be- 
tween both phases are made, can only be achieved in time 
dependent calculations. In two-fluid flow, both phases are 
described by their own continuity and momentum equa- 
tions. Momentum exchange occurs through viscous drag, 
i.e. through gas-grain collisions. The collision rate and the 
momentum exchange per collision depend on the velocity 
of grains relative to the gas. Hence, by flxing the drag 
force, one fixes the relative velocity and the system be- 
comes degenerate. 

In this paper we present our two-fluid time dependent hy- 
drodynamics code. We have selfconsistently included equi- 
librium gas chemistry and grain nuclcation and growth, 
see Section ^ In order not to make assumptions on the 
viscous coupling, we consider, in Section 3.4, the micro- 
physics of gas-grain collisions. Results are given in Section 
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2. Grain drift and momentum coupling 

2.1. Definitions 

The acceleration of dust grains, as a result of radiation 
pressure, leads to an increase in the gas-dust collision rate. 
The viscous drag force (the rate of momentum transfer 
from grains to gas due to these collisions) is proportional 
to the collision rate and to the relative velocity of grains 



with re spcct to the gas. This force is discussed in the next 
section in more detail. The drag force provides a (momen- 
tum) coupling between the gaseous and t he solid phasef]. 



The gas-dust coupling was studied by e.g. Gilman (1972) 



who distinguished two types of coupling. Gas and grains 
are position coupled when the difference in their flow ve- 
locities, the drift velocity, is small compared to the gas 
velocity, i.e. when the grains move slowly through the gas. 
Momentum coupling, on the other hand, requires that the 
momentum acquired by the grains through radiation pres- 
sure is approximately equal to the momentum transferred 
from the grains to the gas by collisions. The situation in 
which both are exac tly equal is cal led full or complete 
momentum coupling. Gilman (1972) stated that, if both 
forces are equal, grains drift at the terminal drift velocity. 
A less confusing term for the same situation was intro- 
duced by Dominik (1992): equilibrium drift. The idea is 
that since the drag force increases with increasing drift 
velocity, an equilibrium value can be found by equating 
the radiative acceleration of the grains and the deceler- 
ation due to momentum transfer to the gas. Note that, 
when calculating the equilibrium value of the drift veloc- 
ity that way, i.e. assuming complete momentum coupling, 
one implicitly assumes that grains are massless. A physi- 
cally correct way to calculate the equilibrium drift velocity 
is to demand gas and grains to have the same acceleration. 

2.2. Single and multi-fluid models 

Various groups have studied the validity of momentum 
coupling, with and without assuming equilibrium drift, 
in stationary and in time dependent calculations. Others 
have just applied a certain degree of momentum coupling 
in model calculations carried out to study other aspects 
of the wind. Wc will give a brief overview of the most im- 
portant of these studies, resulting in the conclusion that 
prior to our attempt, full two-fluid hydrodynamics has 
been presented only twice. Because the meaning of terms 
like "full" and "complete" momentum coupling, "termi- 
nal" and "equilibrium" drift seem to be slightly different 
from author to author, we will first give our own defini- 
tions for three classes of models. 



^ Another momentum coupling is due to the fact that mo- 
mentum is removed from the gas phase when molecules con- 
dense on dust grains. The amount of momentum involved in 
this coupling is also taken into account in our numerical models 
but is many orders of magnitude smaller than the coUisional 
coupling. 



First, single-fluid models are those in which only the mo- 
mentum equation of the gas component is solved. All mo- 
mentum due to radiation pressure on grains is transferred 
fully and instantaneously to the gas. If, e.g., for the calcu- 
lation of grain nucleation and growth rates, a value for the 
flow velocity of the dust component is needed, the dust is 
just assumed to have the same velocity as the gas: drift 



is assumed to be negligible. Hence, in terms of Gilman 



(1972), in single fluid models grains are both position and 



(completely) momentum coupled to the gas. 
The second class is that of the two-fluid models. Here, 
again in terms of Gilman (1972) , grains are not necessar- 
ily position and momentum coupled to the gas. Grains can 
drift at non-equilibrium drift velocities. Hence, grains and 
gas are neither forced to have equal velocity nor forced to 
have equal acceleration. 

The third category of models represents what we will call 
1.5-fluid models. In these models, grains are assumed to 
drift at the equilibrium drift velocity with respect to the 
gas. No assumptions about position coupling are made. 
In other words, gas and grains are equally accelerated but 
do not necessarily have the same velocity. The equilibrium 
drift velocity is calculated by equating the drag force and 



the radiation pressure on the grains, see Dominik (1992) 



or, more accurately, by demanding gas and grains to be 
equally accelerated. Only the momentum equation of the 
gas is solved, the dust velocity is determined by simply 
adding the gas velocity and the equilibrium drift velocity. 

2.3. Stationary models 

Although the above classification for modeling methods 
also applies to stationary models, extra care is needed 
there. When trying to do two-fluid stationary modeling 
one should realize that the condition of stationarity itself 
will also introduce momentum coupling. This can be un- 
derstood as follows. Equilibrium drift is the state in which 
gas and grains are equally accelerated: 

dt dt ^ ' 

The derivative in this equation is a total derivative. Im- 
posing stationarity, the temporal contribution to this total 
derivative vanishes by definition, and Eq.(|l]) reduces to 

The difference between both sides of Eq.(^) can be small, 
especially in the outer layers of the envelope, where the 
velocities reach a more or less constant value. Therefore, 
the occurrence of equilibrium drift in a stationary outflow 
may be partially due to the condition of stationarity itself. 
For this reason, one should be very careful when checking 
the validity of momentum coupling against stationary cal- 
culations. Moreover, in order to make a calculation fully 
self-consistent, no assumptions on momentum coupling 
should be made. Hence, for fully self-consistent modeling, 
time dependent calculations are to be preferred. 
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2.4- Overview of previous modeling 

Examples of single fluid calculations are naturally found in 
studies in which drift and momentum coupling arc not the 
topic of research, e.g. the work of Dorfi & Hofncr (1991) 



and Fleischer et al. (1995). Both perform time dependent 
hydrodynamics, assuming that the influence of drift on 
the aspect of the flow under consideration, dust formation 
and nonlinear effects due to dust opacity, is negligible. 
Th e completeness of momen tum co upling is investigated 
by Berruyer fc Frisch (1983) and by Kriiger et al. (1994) 



The former flrst flnd a (stationary) wind solution under 
the assumption of complete momentum coupling, noticing 
that this assumption causes the two-fluid character to be 
lost. Next, in order to check the validity of their supposi- 
tion, they find a stationary solution for the system, includ- 
ing the grain momentum equation. Both calculations give 
very similar results near the photosphere, from which it is 
concluded that momentum coupling is complete there. Far 
away from the stellar surface (> lOOOi?*), the results are 
different so that momentum coupling is said to be invalid 
there. We too, find that non-equilibrium drift arises far 
away from the photosphere (see Section |4|) . We would like 
to remark, however, that it may not be sufficient to verify 
the validity of complete momentum coupling by compar- 
ing with stationary calculations, see Section e31 



Kriiger et al. (1994) undertook a similar study, which is 
the most realistic stationary two-fiuid calculation up to 
now. It treats the coupled system of hydrodynamics and 
thermodynamics, but also involves chemistry and dust 
formation (simplified by the assumption of instantaneous 
grain formation). Kriiger et al. conclude that momentum 
coupling ca n be assumed to be comple te and therefore dis- 
agree with Berruyer fc Frisch (1983) . We think this may 
be due to the fact that Kriiger et al. run their calculation 
out to about ten stellar radii, whereas Berruyer & Frisch 
compute outwards to several thousand stellar radii. 
According to MacGregor fc Stencel (1992) , who use a sim- 
ple model for grain growth in a stationary, isothermal at- 
mosphere, the assumption of complete momentum cou- 
pling appears to break down for grain sizes smaller than 
about 5 X 10"^ cm. 

Prior to our attempt, time dependent two-fluid hydrody- 



namics was presented by Mastrodemos et al. (1996). They 
conclude that fluctuations on the time scale of the vari- 
ability periods of Miras and LP V (Long Period Variables) , 
200-2000 days, can not persist in the wind. Since they do 
not calculate grain nucleation and growth self-consistently 
but instead assume that grains grow instantaneously and 
have a flxed size, the extreme non-linear coupling between 
shell dynamics, chemistry and radiative transfer (cf. ^edl 
mayr 1997) is not present. Our calculations however indi 



cate that this chemo-dynamical coupling is a main ingre- 
dient to the occurence of variability in the wind. 



similar approach: their models are based on time depen- 
dent, two-fluid radiation hydrodynamics and grains have a 
fixed size. Main emphasis is on the long term variations of 
stellar parameters (L*(t), M(t)), due to the nuclear ther- 
mal pulses, which are included as a time dependent inner 
boundary. It turns out that these large-amplitude vari- 
ability at the inner boundary is not damped in the enve- 
lope and remains visible in the outflow as a pronounced 
shell. 

The calculations presented in this paper aim at combin- 
ing time dependent hydrodynamics with a two-fluid model 
and are suitable for calculating the stellar wind from the 
subsonic photosphere to the supersonic outer layers at 
large distances. We will not take stellar pulsation into ac- 
count because we want to flnd out if the envelope itself 
possesses characteristic time scales. The main goal of this 
work is to get insight in the physical processes underly- 
ing the observed time dependent structures around AGB 
stars. We do not aim at exactly reproducing certain ob- 
servational results and hence will not adjust the stellar 
parameters in order to provide a better flt. 



3. Modeling method 

3.1. Basic equations 

The basic equations for the time dependent description of 
a stellar wind in spherical coordinates and symmetry, are 
the continuity equations. 



dt 



+ ^(''>g^d%d) = Scond,g,d 



and the momentum equations, 
dP 

"dr '^'^^^^'^ JgraVjg "I" VgScondjg 



(3) 



(4) 



dt 



dr 



Steffen and co-workers (Steffen et al. 1997 ; Bteffen et al 
1998; Steffen fc Schonberner 200C ) have a more or less 



frad ~t- ^dragjd fgra,v,d ^g-^condjg (5) 

These equations form a system in which both gas and dust 
are described by their own set of hydro equations (two- 
fluid hydrodynamics). The equations are coupled via the 
source terms. The source term in Eq.(^ represents the 
condensation of dust from the gas, including nucleation 
and growth. Since mass is conserved we have 

-^condjg — -^condjd (6) 

The gas condensation source term is negative due to nu- 
cleation and/or growth of grains. Atoms and molecules 
that condens onto grains take away momentum from the 
gas. This is accounted for in the WgScond.g source terms in 
the momentum equations. 

The momentum equations also couple via the viscous drag 
force of radiatively accelerated dust grains on the gas. 
Since no momentum is lost, we have 

/drag,g /drag,d ('7) 
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The drag force is proportional to the rate of gas-grain 
collisions and the momentum exchange per collision and 
is therefore of the form 

/drag = SdngndTOg|wD|wD (8) 

where Sd is the coUisional cross section of a dust grain 
and vu is the drift velocity of the grains with respect to 
the gas. 

We assume a grey dust opacity and take the extinction 
cross section of the grains equal to the geometrical cross 
section. Then the radiative force is simply 

= (9) 

Radiation pressure on gas molecules is negligible in the 
circumstellar environment of AGB stars. In order to 
determine the temperature structure of the envelope, a 
balance equation for the energy can be added. We do 
not involve the energy structure in the time dependent 
calculation. Also, we do not solve radiation transport. 
Instead, we assume that, throughout the envelope, the 
temperature stratification is determined by radiation 
equilibrium of the gas. This assumption is justified as long 
as the envelope is optically thin to the cooling radiation 
emitted by the dust. The inclusion of an energy equation 
poses no problems, if one wants to spend the computer 
time. 

The model is completed with the equation of state for 
ideal gases. 



3.2. Gas chemistry 

Our hydrocode contains an equilibrium chemistry module 
(iDominik 1992|) which includes H, H2, C, C2, C2H, C2H2 



and CO, and hence is suitable for modeling C stars. 
Oxygen has completely associated with carbon to form 
CO. Due to the high bond energy of the CO molecule 
(11.1 eV), this molecule is the first to form. In absence of 
dissociating UV radiation, CO-formation is irreversible. 
Hence if ec > eo at the time of CO formation, all oxygen 
will be captured in CO and carbon will be available for the 
formation of molecules and dust. Given the total number 
density of H and C atoms in the gas phase, the dissociation 
equilibrium calculation is carried out in each numerical 
time step to give the densities of the molecules mentioned. 
Therefore, bookkeeping of the H and C number densities is 
needed. This requires two additional continuity equations 
of the form of Eq.(|). 

3.3. Grain nucleation and growth 

Once the abundances of the gas molecules are known, 
the nucleation and growth of dust grains can be calcu- 



continuity equations for hydrogen and carbon. The mo- 
ment equations provide the evolution in time of the zeroth 
to third moment of the grain size distribution function. 
Hence, amongst others, the number density and the aver- 
age grain size are known as a function of time. We could, 
in principle, calculate the full grain size spectrum, using 
the moment method, but we limit ourselves to the use of 
average grain sizes. The main advantage of this is that we 
can apply two-fluid, instead of multi-fluid hydrodynam- 
ics, which is obviously computationally cheaper. 



3.4. Viscous gas-grain momentum coupling 

In the absence of grain drift, gas and dust particles will 
collide frequently due to the thermal motion of the gas, 
but no net momentum transfer from one state to the other 
will take place since the collisions are random. If grains 
are radiatively accelerated with respect to the gas, both 
the thermal motion and the acceleration give rise to gas- 
grain encounters, resulting in a net momentum transfer 
from grains to gas. The resulting viscous drag force is de- 



scribed in e.g. Schaaf (1963) 



lated. We use the moment method (Gail et al. 1984; Gail 



fc Sedlmayr 1988), in conservation form (Dorfl & Hofner 



1991 ) . The resulting nucleation and growth rates are used 
to calculate the source terms of Eq. (Q) and the additional 



In the hydrodynamical regime, the time scale on which 
individual gas-grain collisions occur is many orders of 
magnitude smaller than the dynamical time scale. Hence, 
in order to calculate the momentum transfer from grains 
to gas, one needs to sum over many collisional events. 
The strong dependence of the momentum source term on 
the (drift) velocity, via the drag force (Eq.(||)), enables 
rapid changes in the velocities. When applying an explicit 
numerical difference scheme, as we do, it will therefore 
be necessary to take small numerical time steps. Taking 
small, and hence more, time steps involves the risk of los- 
ing accuracy however. In our case, the drag force makes 
the system so stiff that this would lead to unacceptably 
small numerical time steps: a reduction of a factor thou- 
sand or more, compared to the Courant timestep is not 
unusual. To avoid having to take such small steps we per- 
form a kind of subgrid calculation for the drift velocity by 
studying the microdynamics of the gas-grain system. Do- 
ing so, we derive an expression for the temporal evolution 
of the drift velocity during one numerical time step. This 
expression is then used to calculate an accurate value of 
the momentum transfer, i.e. the integrated drag force, in 
one numerical time step. This way, the momentum trans- 
fer rate is determined without making assumptions about 
the value of the drift velocity at the end of the numerical 
time step. Hence, if the momentum transfer is determined 
in this manner a full two-fluid calculation can be done. 
Details of the derivation are given in Appendix 
Another way to go around the problem of course would be 
to assume that the grains always drift at their equilibrium 
drift velocity and to perform a "1.5 fluid" calculation. It 
turns out, however, to be diflicult to determine whether 
or not the assumption of equilibrium drift is justified, c.f. 
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Section 2.3. For a discussion about the comparison of two- 
fluid and "1.5 fluid" calculations see Appendix 



4. Numerical calculations 

4-.1. Numerical method 

The continuity and momentum equations are solved using 
an explicit scheme. A hydrodynamics code was specially 
written for this purpose. It uses centered differencing and a 
two-step, predictor-corrector scheme, applying Flux Cor- 
rected Transport (FCT) ( Boris 1976 ). Second order accu- 
racy is achieved for the single fluid and momentum cou- 
pled ("1.5 fluid") calculations. In the two fluid computa- 
tion we applied, whenever needed. Local Curvature Di- 
minishing (LCD) ( [eke 1991 ), at the risk of introducing 
first order behavior. 



4-. 2. Initial and boundary conditions, grid 

As an initial model for the calculation, a stationary p rofile 
for IRC +10216, kindly provided by J.M. Winters (|Win 



ters et al. 1994 ), was used, see Fig. |^. Stellar parame 
ters of this model are: M, = 0.7Mq, = 2.4 • lO^'L©, 
= 2010K and a carbon to oxygen ratio ec/eo = 1-40. 
The corresponding stellar radius is i?, = 9.20 • lO^'^cm, 
^max = 200i?*. The mass loss rate for the initial model 
is M — 8 ■ lO~^M0yr"^. In order to compare our calcula- 
tions with observations, we extend the computational grid 
to 1287 i?*. Because no initial data is known for the grid 
extension, we simply set the initial values for r > 200i?* 
of all flow variables equal to their value at r = 200i?*. 
As a consequence of this, a transient solution will have to 
move out of the grid before the physically correct solution 
can settle. 

Grid cells are not equally spaced, since a high resolution 
is desirable in the subsonic area but not necessary in the 
outer envelope. The grid cells are distributed according 
to: 



r[n\ — r[n — 1] 



(10) 



r[l] - r[0] 

The number of cells in the grid, rimax; used here is 737 and 
the size ratio q between the innermost and the outermost 
ceh is 318. 

One of the most important aspects of a numerical hydro- 
dynamics calculation is the treatment of the inner bound- 
ary. Since the (long time averaged) mass fluxes throug the 
inner and the outer boundary must be equal, setting the 
inner boundary essentially means fixing the mass loss rate. 
We have, in our calculations, fixed the density and veloc- 
ity in the innermost grid cells, so that the advective mass 
and momentum fluxes (i.e. the first order derivatives of the 
flow variables) through the inner boundary are constant. 
Note that the temperature was constant as a function of 
time as well so that also the pressure will be fixed. In re- 
ality, however, velocity and density will vary with time. 
To account for a variable inflow of mass into the envelope. 



we permit also diffusive inflow of mass. This flow depends 
upon second order derivatives near the inner boundary 
and therefore models quite realistically the cause of matter 
inflow into the envelope. At the inner boundary, the main 
driving term of the wind is not yet active and the velocities 
are very small because newly formed small grains, which 
are very sensitive to radiation pressure, are formed farther 
out. Therefore, the oscillations of the envelope are clearly 
not caused by the implementation of the inner boundary. 
To model the diffusive flux at the inner boundary, we 
could have introduced a separate diffusion term. There 
is no need to do so, however, since our numerical scheme 
involves the calculation of a diffusion term already. This 
diffusion term [numerical viscosity) is part of our finite dif- 
ference scheme and it is locally (i.e. at extrema) required 
to stabilize the centered differencing method. Whenever 
numerical viscosity is not strictly needed to stabilize the 
numerical scheme it will be canceled by an anti-diffusion 
term (Boris 1976). A detailed description of this method 
is beyond the scope of this paper, for details the reader is 



referred to Icke (1991). We want to allow for diffusion at 
the inner boundary. Instead of adding explicitly a diffusion 
term we can simply somewhat reduce the anti-diffusion at 
the inner boundary. That way, not all of the numerical dif- 
fusion is canceled and effectively a diffusive flux is created 
at the inner boundary. 

Although important for the AGB evolution, no stellar pul- 
sations or time dependent luminosities were used. Often, 
in hydrodynamical simulations of late type stars, stel- 
lar pulsations are introduced as a time dependent inner 
boundary condition. In the absence of pulsations, the av- 
erage grain near the inner boundary will be large. Since 
larger grains are less efficiently accelerated by the radia- 
tive force than smaller ones, the stationary inner boundary 
condition will lead to small velocities in the lower enve- 
lope. As a result of the inefficient radiative force on large 
grains, these grains will also tend to drift at high or even 
non-equilibrium drift speeds. To avoid this unwanted be- 
havior, equilibrium drift is imposed in the first 2.8 i?,, also 
in the two-fluid calculation. 

4-3. Calculations 

In order to determine the effect of grain drift on the out- 
ffow, we perform three types of calculation. First, we solve 
the full two-fluid system including gas chemistry, grain 
formation and growth and the continuity and momentum 
equation for both gas and grains. The viscous momentum 
transfer during each numerical time step is calculated by 
integration of /drag over this time step as was presented in 



Section 3.4. Division by the duration of the time step gives 
an expression for /drag that can be inserted in the momen- 
tum equations, Eqs.(|[||). When solving, the left hand side 
of these equations is multiplied by the time step again, so 
that indeed the correct amount of momentum is trans- 
ferred. 
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Fig. 1. Velocity (no drift), gas and dust density, nuclcation rate and average grain radius for the initial profile. 



Next, a 1.5-fluid calculation is performed. Here, the drag 
force is calculated by assuming equilibrium drift in Eq.(^. 
The dust velocity is taken to be the sum of the gas veloc- 
ity and equilibrium drift velocity, according to Eq.(A.34). 
The momentum equation of the dust is not solved. 
Finally, we also perform a single fluid calculation. Here 
too, only the gas momentum equation is solved. The drag 
force exerted on the gas is taken to be equal to the radi- 
ation force on the grains. Now, the velocity of the grains 
is simply set equal to the gas velocity. From the 1.5 and 
single fluid calculations, we expect to learn about the in- 
fluence of (non-equilibrium) drift on the flow, when com- 
paring them to the two fluid calculation. 
All three models were evolved 10^ numerical time steps, 
which amounts to 9.71 • 10^°, 1.67 • 10" or 3.14 • 10" sec- 
onds, depending on the model. 



4-4- Results 

Fig. 1^ shows the mass loss rate at i? = 100, 500 and 1000 
i?, as a function of time for the three calculations. The 
first 150 years of output in the 500 i?* plot and the first 
800 years in the 1000 i?* plot show the passing of the 
transient solution. This is a result of extending the grid 
from 200 i?* in the initial profile to 1287 i?, in the calcu- 
lation, the flow needs some time to reach the additional 
gridpoints. 

Both the 1.5 and the two-fluid model show quasi-periodic 
oscillations. From plots which cover a longer time inter- 
val (not shown here) we infer that the variations in the 
mass loss rate in the single fluid calculation behave quasi- 
periodically as well, on a time scale of a few thousand 
years. An immediate conclusion from this is, that the pres- 
ence of grain drift is important for variations of the mass 
loss rate. 

The time between two peaks in the mass loss is approxi- 
mately 200 to 350 years for the 1.5-fluid model, and about 
400 years for the two-fluid model. Both numbers lie nicely 
in the range of the separation of 200-800 years between 
the shells that Mauron & Huggins (1999)1 observed in IRC 
-M0216. 

In all three calculations we see that the short time vari- 
ations that are pres ent at 100 R^, have disapp eared far 
away from the star. Mauron fc Huggins (2000l note that 



this "wide range of shell spacing, corresponding to time 
scales as short as 40 yr (close to the star) and as long as 
800 yr", should be accounted for in a consistent model. 
This poses no problems, since the disappearance of the 
smaller scale structures is simply due to dispersion and 
hence will appear in any flow in which perturbations do 
not propagate with exactly the same speed. 
The fact that the two-fluid calculation shows less varia- 
tions on short times scales than the 1.5-fluid model may 
be due to the more first order character of the former (as a 
result of the LCD term, see Section 4.1). We shall see that 
in the two-fluid calculation, in large parts of the envelope, 
grains move at their equilibrium drift velocity. The time 
averaged mass loss rate, estimated from Fig. lies around 
M — 1 ■ lO^^'Moyr^^. The fact that this is somewhat 
higher than the mass loss rate of the initial model indi- 
cates that indeed the diffusive flux at the inner boundary 
has contributed, see Section |4.2| . Our limited implemen- 
tation of the radiative force (we use a grey dust opacity 
and take the extinction cross section of the grains equal 
to the geometrical cross section) causes the velocities in 
our calculation to be higher than the velocities in the ini- 
tial model. Using a lower value for the stellar luminosity 
(e.g. using the core mass-luminosity relation) has proven 
to immediately lower the outflow velocity and hence the 
mass loss rate. 

Figs. H and ^ show, for the 1.5 and the two-fluid model, 
the gas and dust velocities and densities, as a function of 
radius and time. Throughout the whole grid, the fluctua- 
tions occurring in the two-fluid calculation are more reg- 
ular that those in the 1.5-fluid model. The velocities of 
gas and dust in the momentum coupled calculation reach 
values that are up to 25% higher than in the two-fluid cal- 
culation. In the latter, matter is less accelerated than in 
the former, especially for radii larger than about 2 • 10^^ 
cm. Probably, this is a result of non-equilibrium drift, 
which starts to appear around this radius (see Fig. ^. 
Non-equilibrium drift occurs when the time needed by a 
grain to reach its equilibrium drift velocity is long com- 
pared to the dynamical time scale. During a period of 
non-equilibrium drift, the gas is not being maximally ac- 
celerated and both gas and dust velocities will be lower 
than in a phase of equilibrium drift. 

The gas density structure (Fig. ^) for the 1.5-fluid and the 
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Fig. 2. From top to bottom: Mass loss rates for single fluid (no drift, gas and grain have equal velocity, "position 
coupling"), 1.5-fluid (equilibrium drift, gas and grains have equal acceleration, "momentum coupling") and two-fluid 
(no assumptions on drift, no coupling imposed) calculations, for R = 100, 500 and 1000 i?*. Note that the first 150 
years of output in the 500 i?* plot and the first 800 years in the 1000 i?* plot show the passing of the transient solution 
due to the extension for the calculational grid w.r.t. the intial model. 



two fluid calculation look similar. The main difference is 
that short time scale variations are present in the lower 
regions of the former, whereas large scale effects domi- 
nate the latter. The density structure plots for the dust 
show another difference: the perturbations in the 1.5-fluid 
flow appear as local increments of the density but in the 



two component flow the variations rather look like dips 
in the average proflle. Maximum outflow density for gas 
and grains are in phase in the two-fluid model though, 
the "dust pulse" is significantly broader than but centered 
around the maximum in the gas outflow. This is not just 
the case in the upper parts of the envelope, where non- 
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Fig. 3. Gas and dust velocities as a function of radius and time for the 1.5 and the two-fluid model. 



equilibrium drift is present, but also for smaller radii. 



mined in the subsonic region (see e.g. Lamers & Cassinell 



In Figs . 5 and 3 we plot a series of snapshots, displaying 
the evolution of various flow variables during one insta- 
bility cycle for the 1.5 and the two-fluid model. For the 
1.5-fluid calculation the drift velocity is, by deflnition, al- 
ways equal to its equilibrium value, which shows a time 
dependent behavior. In the two-fluid flow we find that the 
drift velocity, out to approximately 10^^ cm, equals the 
equilibrium value. At larger radii, small deviations from 
equilibrium drift are detected. 

We want to stress that the fact that we see equilibrium 
drift in the lower and intermediate regions of the two com- 
ponent model only implies that equilibrium drift is estab- 
lished on a time scale shorter than the dynamical time 
scale. It does not however exclude the possibility that 
non-equilibrium drift occurs on shorter time scales, see 
Appendix 

4-. 5. The origin of the mass loss variability 

To investigate what causes the variability we will step 
through the frames of Fig. ^ for the two-fluid calculation. 
Thereafter, we will discuss the differences with the 1.5 
fluid model. The mass loss rate of a stellar wind is deter- 



1999), therefore in the following, when investigating the 



mechanism underlying the variability, we focus on this re- 
gion, unless explicitly mentioned. 

In Fig. H, first frame, we see that the onset of the mass 
loss variability is the situation in which the dust has a 
velocity that is significantly higher than the gas velocity. 
This means that the residence time of a grain in the parts 
of the envelope where grains can grow is relatively short 
so that the average grain size will be on the small side. 
The smaller the grain, the more efficient radiation pres- 
sure will be, since small grains have a large surface to mass 
ratio and since we have assumed that the grain extinction 
cross section equals the geometrical cross section. Hence, 
radiative acceleration of grains is efficient and the veloc- 
ity of the small grains increases further. Because position 
coupling is not imposed, the gas velocity can stay low 
and the drift velocity increases. Meanwhile (frames 2 and 
3), the average grain radius decreases, grain acceleration 
becomes more efficient, the dust velocity grows, grains be- 
come smaller, and so forth. Also, the total mass density 
of the dust component in the innermost region decreases. 
When the grain radius in the subsonic region drops below 
a certain critical value, momentum transfer from grains to 
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Fig. 4. Gas and dust densities as a function of radius and time for the 1.5 and the two-fluid model. 



gas becomes efficient and the gas is accelerated (frame 4) . 
This results in an increase of the gas density and hence 
of the number density of condensible particles. Since the 
grain nuclcation rate is extremely sensitive to the molecu- 
lar abundances, this results in an immediate increment of 
the nucleation rate (frame 4). The new production of con- 
densation kernels leads to a further decrease of the average 
grain radius and an increase of the total grain mass den- 
sity. Due to the large abundance of small grains, radiative 
acceleration and the transfer of momentum from grains to 
gas are very efficient, so that both gas and grains move 
out with high velocities (frames 5-8). On their way out, 
the small grains concentrate in a narrowing shell, since the 
decrease of the average grain radius in time coincides with 
an increase of their velocity. The gas develops a shell at 
the same time, as a result of the forming shock. The nor- 
mal, Parker-type, stellar wind profile is now visible. We 
will refer to this phase as the "fast phase" (frame 5-9). 
Though not very clear from the figure, at the same time, 
a rarefaction wave moves in the opposite direction, lead- 
ing to a decrease of the gas density, and of the number 
densities of the condensible species, below the sonic point. 
Although the density decrease is not so big, the nucle- 
ation rate reacts instantaneously (frames 9-13), showing 



a strong decrease traveling from the sonic point inwards. 
Hence, the passing of the rarefaction wave is immediately 
visible in the increase of the average grain radius because 
the production rate of new small grains decreases (frames 
9-13). This illustrates the enormous sensitivity of the nu- 
cleation rate on the densities. The gradual increase, in 
time, of the average grain radius, brings about a less effi- 
cient radiative acceleration of the dust, hence a decrease 
of the grain velocity and a further increase of the grain 
radius, and so forth. This wc will call the "slow phase" of 
the variability cycle (frames 10-14 and 1-4). Due to the 
larger grain size, the momentum transfer between grains 
and gas becomes less efficient, resulting in larger drift and 
dust velocities (frame 14). This brings us back to the sit- 
uation in the first frame. 

Crucial in the process of shell formation as described 
above are the two "turn-around" points, at which the nu- 
cleation rate starts to increase and decrease. First, at the 
end of the fast phase, the passage of the rarefaction wave 
triggers the end of a period of high nucleation rate. In 
the slow phase the gas-grain coupling has becomes less 
efficient, due to the larger average grain size. Grains then 
reach a higher drift velocity, become smaller and will again 
transfer their momentum efficiently to the gas, so that 
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Fig. 5. 1.5-fluid model. First column: gas and dust velocity (dashed line). The dot denotes the location of the critical 
point. Second column: gas and dust density (dashed line). Third column: drift velocity. Fourth column: average grain 
radius and grain nucleation rate (dashed line). The frames show (from top to bottom) the flo'w profile at 0, 30, 56, 81, 
105, 132, 164, 197, 225, 252, 280, 310, 352 and 404 years after the first frame. 
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Fig. 6. Two-fluid model. First column: gas and dust velocity (dashed line). The dot denotes the location of the critical 
point. Second column: gas and dust density (dashed line). Third column: drift velocity (dashed line) and equilibrium 

drift velocity (full line). Fourth column: average grain radius and grain nucleation rate (dashed line). The frames show 
(from top to bottom) the flow profile at 0, 32, 63, 93, 116, 131, 144, 156, 170, 189, 214, 245, 284 and 330 years after 
the first frame. 
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the latter can accelerate, increasing the density. This gives 
rise to favourable circumstances for grain nucleation again. 
Clearly, the behavior of the system during the slow phase 
is dominated by the existence of grain drift. This imme- 



stifF and reacts violently to all kinds of changes, we think 
that it is rather unlikely that the observed circumstellar 
shells are indeed complete. It is intriguing to see that this 



idea is supported by the recent observations by Mauron 



diately explains why variability in the mass loss rate in a 
single fluid system is less well regulated (see Fig. |^). 
When comparing Fig. ^ and Fig. |[ the absence of the slow 
phase in the variability cycle in the latter strikes the eye. 
This can be attributed to the imposed equilibrium drift in 
the 1.5-fluid flow. In the two-fluid system the drift veloc- 
ity is directly influenced by the dynamics. In the 1.5-fluid 
model, however, the (equilibrium) drift velocity is only 
indirectly determined by the dynamics, namely via the 
(number) densities and the grain size. The fact that the 
variability character is still observed in this calculation is 
a consequence of the fact that the drift velocity, although 
not actively, does change as a function of time, in combi- 
nation with the extreme sensitivity of the nucleation rate 
to the density and of the dynamics, via the drag force, 
on the grain size, and density. The sensitivity of the sys- 
tem is well visible in Fig. any variation of the densities, 
grain size and nucleation rate is hardly visible (also be- 
cause they are plotted logarithmically, ranging over many 
orders of magnitude) but the resulting variations in the 
velocity field are clearly present. 

4-6. Comparison with observations 

To enable a qualitative comparison of our results with 
recent observations of IRC -1-10216 (Mauron & Huggins 
1999; 2000), we have produced Fig. [?[ T he left frame is 
adapted from Mauron & Huggins (1999) (their Fig. 3). It 
shows the composite B + V image of IRC -1-10216, with 
an average radial profile subtracted to enhance the con- 
trast. We compare this image with the dust column den- 
sity as a function of radius for a number of snapshots in 
our calculation. The size of our computational grid (ex- 
tended to 1287 i?*) corresponds with the field of view of 
the observational image (131" x 131") and a distance of 
120 pc. We too, have subtracted an average radial density 
profile to enhance the contrast. Comparing dust column 
density to the observed intensity makes sense, since in the 
optically thin limit, the observed intensity, due to illumi- 
nation by the interstellar radiation field is proportional 



Huggins (2000), which show that most shells, although 



to the column density along any line of sight (Mauron & 



Huggins 2000). We used the results of the 1.5-fluid com- 
putation to produce Fig. ^ because there the short time 
scale structures are visible, whereas they are suppressed in 
the two-fluid model because the latter isn't always second 
order accurate. Note that the fact that in our calculated 
images all shells appear to be perfectly round is simply 
due to our assumption of spherical symmetry. The two 
dimensional plots were produced by simply rotating the 
spherical symmetric profile. In view of the fact that our 
calculations indicate that the chemical-dynamic system 
that regulates the behavior of the envelope is extremely 



they may extend over much larger angles at lower levels, 
are prominent over about 45°. 

As was mentioned before. Fig. ^ only offers a qualitative 
comparison with the observations. It can, however be used 
to establish that the spacing of the shells, small scale struc- 
ture inside, large scale structure outside, is similar in the 
observations and calculations. This, is not surprising how- 
ever, since merging of shells of various widths is due to 
dispersion, as was mentioned in Section 4.4. 



4-.7. The timescale of mass loss variations 

The characteristic time scale of the variability corresponds 
to the time needed by the rarefaction wave to cross the 
region between the sonic point and the innermost point of 
the nucleation zone. The width of this region is, depending 
on the phase, a few times lO^'' to 10^^ cm. The velocity 
of the rarefaction wave equals the gas velocity minus the 
local sound velocity and is typically a few times lO'^ to 
10^ cm s~^, also depending on the phase of the variability. 
The resulting time scale is roughly 50 to 500 years, which 
indeed corresponds to the time separation between two 
maxima in the mass loss rate in our calculation. 



4-.8. Discussion 

We found that the fact that the average grain size reacts 
strongly to the density structure is an essential ingredi- 
ent for the fo rmation of variability in the out flow. This 
explains why Mastrodemos et al. (1996) and Steffen fc 
^chonbcrncr (2000)| , who also performed time dependent, 
two-fluid computations, but did not take into account self 
consistent grain growth, did not encounter mass loss vari- 
ations in the outflow. 

Also, grain drift occurs to be essential for variations in the 
mass loss rate. If grains can drift with respect to the gas, 
they can form regions of higher (or lower) density and/or 
size independently from the gas. 

Periodic variability in the mass loss rate occurs in both the 
1.5-fluid and the two-fluid calculations, because grains 
are allowed to drift in both cases. Both calculations give 
somewhat different results, though. Probably, assuming 
equilibrium drift a priori, as was done in the 1.5-fluid 
computation, influences the results, even if the grains in 
the two-fluid model turn out to drift at the equilibrium 
drift velocity as well. There are two reasons for this. First, 
the fact that equilibrium drift has established itself at the 
end of a numerical time step, does not mean that there 
has been equilibrium always during this specific time step. 
Hence, integration of the drag force over the time step pro- 
vides a better value of the momentum transfer than mul- 
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Fig. 7. Upper left frame: Composite B + V image of IRC +10216, with an average radial profile subtracted to enhance 
the contrast (adapted from Mauron & Huggins (1999)). Note that a few patches in the image are residuals of the 
removal of the brightest background objects, and these should be ignored. Other frames: series of snapshots of our 
1.5-fluid calculation. Plotted is the dust column density, also with an average radial profile subtracted. The average 
radial profiles are calculated for each snaphot separately, hence the slight difference in color from plot to plot. The 
theoretical profiles are shown for ages 44, 118 and 211 years with respect to the first frame in Fig. ||. The size of our 
computational grid corresponds with the field of view of the observational image (131" x 131") and a distance of 120 
pc. 



tiplication of the drag force with the duration of the time 
step, c.f. Appendix ^. Second, the value of the equilibrium 
drift velocity in the 1.5-fluid calculation is indirectly de- 
termined by the dynamics, whereas in the two-fluid case 
there is a direct influence. Also, the fact that the 1.5-fluid 
calculation is second order accurate, but in the two-fluid 
calculation this level of accuracy is not always achieved, 
will lead to differences in the results. 
We have not taken into account radiative transfer to solve 
the energy structure in the envelope. Also, we used a grey 
absorption coefficient in the radiative force and we did 
not calculate the grain temperature. These are severe lim- 
itations of the model. However, we believe that they do 
not influence the general conclusion that dynamics and 
chemistry together can lead to time dependent structures. 
It is more likely that taking into account the tempera- 
ture structure determined by the optical properties of the 
grain population will make the variability even more pro- 
nounced. This is inferred from previous calculations by 
Fleischer et al. (1992)| in which the interaction between at- 
mospheric dynamics and radiative transfer was solved, im- 



posing a time dependent inner boundary. Recently, Win- 
ters et ^1. (2000) performed similar calculations, also with- 



out the piston at the inner boundary. Their results also in- 
dicate that the coupling between the sensitive grain chem- 
istry and the dynamics can lead to variability in the wind. 
The role of the inner boundary in calculations as presented 
here is extremely important. It is possible to generate wind 
variability using a time dependent inner boundary. We did 
not do this: the inner boundary that we have used was cre- 
ated to have as little influence on the results as possible. It 
consists of a fixed advective flux which can be modified by 
a diffusion term. The diffusive contribution to the flux is 
proportional to the gradients of the flow variables near the 
inner boundary, i.e. it is not externally prescribed. This is 



a realistic approach, since the inner boundary is located 
in the subsonic regime, where communication with lower 
layers is still possible. In this respect a completely fixed 
inner boundary would be less realistic. 
We have referred to the quasi-periodic structure in our 
models as "shells". In order to prove that the structure 
is truly created in the form of spherical shells one should 
perform three dimensional hydrodynamics. Higher dimen- 
sionality will be a topic of future research. Shell structure 
is observed around only a small number of Post- AGB ob- 
jects and PNe. It is possible that the majority of objects 
doesn't have shells. A stationary wind can definitely exist 
if for some reason the equilibrium drift velocity is rela- 
tively low. This can be the case if the luminosity of the 
star is low. This will limit the mutual motion of both flu- 
ids and hence the value of the gas to dust density ratio so 
that the outflow will remain more smooth. 

5. Conclusion 

Our calculations suggest that the sensitive interplay of 
grain nucleation and dynamics, in particular grain drift, 
leads to quasi-periodic winds on the AGB. The charac- 
teristic time scale for the variability corresponds to the 
crossing of the subsonic nucleation zone by the rarefaction 
wave. This time scale also matches recent observations of 
IRC -hl0216. 

More generally, we would like to stress that two-fluid 
hydrodynamics is important in order to reach self- 
consistency of the modeling method since the validity of 
the assumption of equilibrium drift is hard to check. // 
equilibrium drift is applied, it should be calculated by de- 
manding the grains and the gas to be equally accelerated, 
rather than by equating the drag force and the radiation 
pressure on grains, because grains do have mass. 
Observations also imply that gas and grains may not be 
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Fig. A.l. Evolution towards equilibrium drift within a 
single time step: because it may take some time to estab- 
lish equilibrium, calculating the momentum transfer by 
simply multiplying the drag force (which is proportional 
to Uq) with At overestimates the momentum transfer. The 
difference between the exact calculation and the equilib- 
rium calculation increases with the time required to estab- 
lish equilibrium drift and is represented by the dark color 
in the figures. 



spatially coupled (Sylvester et al. 1999) and that varia- 
tions in the gas to dust ratio in the outflow may arise 
( pmont et al. 1999| ). 
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Appendix A: Calculating the drag force 

To derive an expression for the drag force, we need to 
know about the time evolution of the drift velocity. The 
gas-grain system will always evolve towards a state in 
which grains drift at the equilibrium drift speed, hence 
in which gas and grains undergo the same acceleration. 
If, or how rapidly this state is reached depends on the 
time needed to establish the equilibrium relative to the dy- 
namical time scale. If one assumes that equilibrium drift 
is always valid, the momentum transfer in a numerical 
time step can simply be calculated by using the equilib- 
rium value of the drift velocity in Eq.(^) and multiplying 
the drag force by the duration of the time step. How- 
ever, if, during a fraction of the numerical time step, the 
drift velocity is lower than the equilibrium value, assum- 
ing equilibrium drift when calculating the drag force will 
overestimate the momentum transfer. This is illustrated 



in Fig. A.l. Although the error for a single time step may 
be very small, the implications may be large for the time 
dependent calculation. Note that, when assuming equilib- 
rium drift, one fixes the value of the drift velocity so that 
the gas and the dust velocities are no longer independent 
flow variables. Therefore, when calculating the momen- 
tum transfer assuming equilibrium drift one is forced to 
do a 1.5-fluid calculation rather than a full two-fluid cal- 
culation. We will, hereafter, derive an expression for the 
time evolution of the drift velocity. With this expression 
we can calculate the momentum transfer as the integral of 
the drag force over the numerical time step. No assump- 




Fig. A. 2. Evolution towards equilibrium drift for various 
initial drift velocities. Upper panels: go, tot > wd > 0, 
lower panels: 5D,tot < — > < 



tions about the final drift velocity need to be made and 
the derived expression can be used in a full two-fluid cal- 
culation. 

It is important to note that even if we find equilibrium 
drift in the two component calculation this does not im- 
ply that it would have been justified to assume equilibrium 
drift a priori. This can be seen from Fig. A.l. In both the 
first and the second panel equilibrium drift is established 
within the duration of the numerical time step, At, i.e., 
in both cases the output of the hydrodynamics indicates 
equilibrium drift. Assuming equilibrium drift throughout 
the time step would however only slightly overestimate 
the momentum transfer in the first panel whereas is the 
second panel the difference between the exact integral of 
the drag force over the time step and equilibrium approx- 
imation would be much bigger. 



A.l. An analytical expression for the momentum transfer 
rate 

In this section we will derive an expression for the time 
evolution of the drift velocity. Using this expression we 
can calculate the rate at which momentum is transfered 
from grains to gas. 

Fig. A.S shows the six possible cases for reaching equilib- 
rium drift. Note that both the initial drift and the equi- 
librium value can be negative if the grains are less accel- 
erated than the gas. We assume that the gas-grain inter- 
actions are completely inelastic. Furthermore, we assume 
that after a collision with a grain, a gas particle shares the 
acquired momentum with the surrounding gas instanta- 
neously (thermalization) . This is realistic, since the mean 
free path of gas-gas collisions is very small compared to 
the mean free path for gas-grain encounters. We will not 
take into account thermal motion because this enables us 
to derive an analytic expression for the drag force. This 
will result in a somewhat lower momentum transfer in 
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the subsonic region. Farther out, the drift velocity of the 
grains will dominate the collision rate anyway. 
First, consider the motion of an individual gas particle 
between two subsequent collisions with a grain: 

^ Wg + ge,totSt + (A.l) 

ng mg 

Here, Vg is the velocity of the particle, after the previ- 
ous collision, gg,tot is the total acceleration due to gravity 
and the pressure gradient (but not the drag force), St is 
the time interval between two collisions. The last term 
represents the increase in the velocity as a result of the 
encounter with the grain, and the (instantaneous) redis- 
tribution of the momentum amongst the gas. Ap is the 
amount of momentum transferred in a single gas-grain 

collision, 

mgmd 



Ap = 



(A.2) 



mg + md 

where ud is the velocity of a grain with respect to the 
gas immediately before the collision, rrig^d are the masses 
of a gas particle (i.e. the mean molecular weight) and the 
(average) grain mass. A similar equation for the dust grain 
is 

Vd+ gd,totSt - — (A. 3) 

md 

The drift velocity after a collision, vd, can now be ex- 
pressed in terms of the drift velocity immediately before 
the encounter, uu, as follows: 

vb = fluY) (A. 4) 

in which 

n = P^"^-" - P'^^l (A.5) 

Pg(mg -Hmd) 

wd = Wd - Ug = Wd - t^g + (5d,tot - 5g,tot)^i (A. 6) 

In the following, we will write gD.tot for the relative ac- 
celeration, (7d,tot — 5g,tot- The "mean free travel time", (5t, 
of a grain can be found by solving the quadratic equation 
for the mean free path. A, of a grain 

A = v^St + igo.tof^i^ 

Note that the mean free path can become negative if the 
initial drift velocity, wd, and/or the relative acceleration 
ffD.tot is negative. If grains are not significantly accelerated 
between two subsequ ent collisions with gas particles, i.e. 
if ud 3> go.tot^i, Eq.(A/?) simply becomes 
A = vxySt (A.8) 
so that St = A/ud. On the other hand, if the acceleration 
of a grain between two collisions is so large that its initial 
(drift) velocity is negligible, Eq.(A.7) reads 

A = ^ffD,tot<5i' (A.9) 



and St = A/2A/(7D,tot- The boundary between the two 
regimes lies at the drift velocity for which 2wd = 5D,tot^i- 



With St given by the solution of Eq.(A.7) we find that if 

\M<\\fM^t (A.IO) 
Eq. 



3|) can be used instead of Eq.( A.7 ). In the current 
context of dust forming stellar winds, the quantity will 



always be nearly equal to unityQ, so that vy^^ \ \/AgD,tot- 
Hence, the zone in velocity space where grain acceleration 
is significant is extremely narrow. If the drift velocity is 
zero at some time (see e.g. Fig. [A.2| .c,f), it follows from 
Eq.(^!|), ( |a!6| ) and (|Al9| ) that the drift velocity will be 
larger than 2\/A5D,tot after a single collision unless < 
l/VS- This implies that we can safely apply Eq.( |A.8| ) for 
all values of • 

In the following we will present a method to derive an 
expression for the momentum transfer, which applies to all 
possible scenarios (see Fig. A.2) to reach equilibrium drift. 
We limit ourselves to the derivation for the case ^D.tot > 
(Fig. A.2.a,b,c), the derivation for negative acceleration is 
analogous. 



Application of Eq.(A.8) and Eq.(A.6) in Eq.(A.4) gives 



rise directly to a recurrence relation for vy^: 

vu{t^+i) = r!uD(i»+i) = n (vo{ti) + ^^^] (A.ll) 

From this, and St given by Eq.( |A.8| ), a differential equation 
for the drift velocity as a function of time can be derived: 

+ (A.12) 
This equation can be easily solved for i(wD), 
t(wD = , 

^n{n - i)gx 

arctan ( i^^l^L^ _ 



. ( (""f> D(o) 

arctan — , 



(A.13) 



where g stands for gD,tot- 

First, consider the case where wd(0) > (and g > 0). 
In this case the mean free path A will always be positive. 
Because fl is always smaller than unity and A and g have 
equal signs this is rewritten as 
A 



(A.7) tivu) = 



arctanh 



n)g\ 

( (i-n)vT,{t) ^ 
\^n{i-n)g\^ 

' (i-r!)wD(o) ' 



arctanh | -y-i^y^ ^ . (k.X^) 

This expression can be simplified by realizing that from 
Eq.(A_.ll) it follows that the equilibrium drift velocity is 
given by 



-M 





and that the equilibration time scale is 
1 



' cq 



(A.15) 



(A.16) 



Vr!(l - ^)gl\ 

E.g. for a typical dust to gas mass ratio pd/pg = 1.0 x 10~^ 
and for grains consisting of 10^" momomers (md/mg — 1.0 x 
10^°) we find f« 1 - 10"^° 
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so that 

tivu) = 



arctanh 



arctanh 



arctanh 



vujt) 



wd(0))ud 



(A.17) 



(A.18) 



Note that addition of the arctanh terms causes the ex- 
pression to be vahd for initial values wd(0) > vb (see 
Fig. A.2.b) as well. Inversion leads to an expression for 
the drift velocity as a function of time: 

VD(0) + WDe(t) 



(A.19) 



with 

Q{t) = tanh(i/Tcq) (A.20) 
The drag force (density) is the product of the number 
of gas-grain collisions per unit volume and time and the 
momentum transfer per collision. In Eq.(^), the amount 
of momentum transfer in a single collision was simply as- 
sumed to be TOgWD, now we use the more accurate form 
for Ap which follows from Eqs.(^), ( |A^ ), With 
A = l/Sd^g we then find 

/drag = SdPg IwdI^'D (A. 21) 

rig — rid 

The standard way to calculate the amount of momentum 
transfer per numerical time step is simply multiplying the 
drag force with the duration of the time step. Now that we 
have derived an expression for the drift velocity as a func- 
tion of time we can calculate the momentum transfer more 
accurate, by integrating Eg. ( A. 21 ), assuming ng_d,TOg,d 
are constant: 

rigrid 



/dragdi — SdPg 



rid 



wd(0) tanh(T/Tcq) 



f vuiO) 

\ vu UD(0)y V"d(0) tanh(T/rcq) -I- WD 



(A.22) 



If the initial drift velocity and the total acceleration have 
opposite sign (wd(0) < 0, .g > 0, see Fig. A^.c) the integral 
representing the total momentum transfer is split into two 
parts, 

(.t(l>D=0) nT 

/dragdt = / /dragdt + / /dragdt (A. 23) 

Jo J t(vD= 0) 

where t{vB = 0) follows from Eq.( A.l^ ): 
tivB = 0) 



y/n{n - i)g\ 

arctan ( S^l^^^] (A.24) 

\^ - i)g\ J 

Note that the mean free path of a grain. A, is negative as 

long as the drift velocity is negative. The second term in 

Eq.( A.23| ) is calculated as in the case i;d(0) > 0, simply 

taking wd(0) = 0. In order to compute the first term, 

Eq.( [A.13| ) is inverted. We find 

... _ vu{0)+VB&{t) 

VD{t) = VB- (A.25) 

VB - VB(0)Q'[t) 



in which 

e'{t) = tan(Vr' ) 



Vb 



■cq 



n 



n- 1 



(A.26) 
(A.27) 
(A.28) 



Inserting this into Eq.( A.21 ) and integrating over the in- 
terval t — 0, t{vD = 0), we obtain 

^t(l>D=0) 



/dragdt — — SdPg 



ngnd 



Tig - ?ld 

vb{0) 



'cq 



arctan 



VBiO) 



(A.29) 



Note that the minus sign accounts for the fact that the mo- 
mentum transfer contains an integral over Iwd^d rather 
than an integral over . Finally, for the complete integral, 
Eq.( |A.23D , we find 

TlgTld 



/dragdt — SdPg 



rid 



-'''cqVB 



tanh 



' cq 



arctan 



' cq 



^d(O) 
Vb 



VBiO) 



VB 



(A.30) 

As was to be expected Eq.( A.2 j ) and Eq.( A.3(i| ) arc equal 
if 1^0(0) = 0. 

Similar expressions for the total momentum transfer can 
be calculated in the case of negative total acceleration (see 
Fig. ^d,e,f). 

The above formulations for the momentum transfer, in 
which no assumptions about the value of the drift veloc- 
ity or the completeness of momentum coupling have been 
made, can be used as source terms in the momentum equa- 
tions. 

A. 2. Calculation of the equilibrium drift velocity 

We have used the terms equilibrium drift velocity and lim- 
iting velocity as equivalent. Here, we will show that both 
are indeed the same. We equate the acceleration of the 
gas and the dust, rather than equating the drag force and 
the radiation pressure of grains. In the latter case one im- 
plicitly assumes that grains do not have mass whereas the 
former leads to a general expression for the equilibrium 
drift velocity. 

From the equation of motion of a gas element, 

_ „ I /drag (\ '^^\ 

^-ffg,tot + ^ (A.31) 
and its counterpart for a grain, 

_ /drag , . „„n 

— r- — .9d,tot (A.6Z) 

at pd 

we find that grains and gas are equally accelerated, and 
hence the drift velocity has reached its equilibrium value, 
if 



5D,tot — 



Pd + Pg 
Pd£g. 



drag 



With Eq.( A.21 ), the equilibrium drift velocity is 



Vb = 



I mdjug - U d) 
Pd + Pg 



-5D,tot 



A 



(A.33) 



(A.34) 
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Thus, we have now derived an expression for the equi- 
librium drift velocity without having to assume complete 
momentum coupling. This expression is indeed the same 



as Eq.(A.15), which represents the limiting drift velocity. 
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